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> ■ ABSTRACT 

\^ , We study the influence of the matter content of extragalactic jets on their mor- 

phology, dynamics and emission properties. For this purpose we consider jets 
of extremely different compositions including pure leptonic and baryonic plas- 
' mas. Our work is based on two-dimensional relativistic hydrodynamic simula- 

tions of the long-term evolution of powerful extragalactic jets propagating into 
a homogeneous environment. The equation of state used in the simulations 
accounts for an arbitrary mixture of electrons, protons and electron-positron 



C/3 

^ ' pairs. Using the hydrodynamic models we have also computed synthetic radio 



maps and the thermal Brcmsstrahlung X-ray emission from their cavities. 

Although there is a difference of about three orders of magnitude in the 
temperatures of the cavities inflated by the simulated jets, we flnd that both 
the morphology and the dynamic behaviour are almost independent on the 
assumed composition of the jets. Their evolution proceeds in two distinct 
epochs. During the flrst one multidimensional effects are unimportant and 
the jets propagate ballistically. The second epoch starts when the first larger 
vortices are produced near the jet head causing the beam cross section to 
increase and the jet to decelerate. The evolution of the cocoon and cavity is in 
agreement with a simple theoretical model. The beam velocities are relativistic 
(r ~ 4) at kiloparsec scales supporting the idea that the X-ray emission of 
several extragalactic jets may be due to relativistically boosted CMB photons. 
The radio emission of all models is dominated by the contribution of the hot 
spots. All models exhibit a depression in the X-rays surface brightness of the 
cavity interior in agreement with recent observations. 
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1 INTRODUCTION 

The standard model for powerful jets associated with extragalactic radio sources (Blandford 
& Rees 1974) assumes that the energy radiated in the radio lobes of such sources is produced 
in the active nuclei of their host galaxies, the central engine being a supermassive black hole 
surrounded by an accretion disc. The accretion process fuels a couple of twin supersonic 
jets which transport away bulk kinetic energy from the neighbourhood of the central black 
hole to the lobes, i.e., from scales of the Schwarzschild radius of the supermassive black hole 
Rs = 2GMbh/c^ = 3 ■ lO"'^''^ ■ (Afbh/IO^MQ) cm, to kiloparsec scales. This bulk kinetic energy 
is dissipated by shocks within the beam and (mostly) at the jet terminal shocks where 
electrons are accelerated and radiate via synchrotron and inverse Compton mechanisms. 
Typical lifetimes and kinetic powers of powerful radio sources are 10^ y and 10^^ — 10^^ erg 
s~^, respectively (Rawhngs & Sunders 1991; Daly 1995). 

Among the problems that still remain open even after more than 30 years of research is 
the composition of extragalactic jets. Within the standard model (Blandford & Rees 1974) 
jets are made of a relativistic plasma that contains rclativistic electrons and thermal protons 
(ep-plasma). On the other side, based on equipartition arguments Kundt & Gopal-Krishna 
(1980) claim that extragalactic jets consist of beams of extremely relativistic electrons and 
positrons (e^-plasma) with almost no ions. Ultimately, the composition of jets is tightly 
related to their formation mechanisms. As discussed by Celotti & Blandford (2001), electro- 
magnetically dominated outflows, as those generated by the extraction of spin energy of the 
black hole, will become pair dominated jets with a low baryonic pollution. Jets generated 
from the accretion disk by hydromagnetic winds will be made of baryonic plasma. Probably, 
both processes are operating simultaneously in nature. Sol, Pelletier & Asseo (1989) propose 
a two-flow model where the jet consists of a beam of relativistic particles (pair plasma) sur- 
rounded by a Newtonian or mildly relativistic (ep-plasma) wind accelerated from the disk. 
It is also possible that ep-jets become pair-loaded later on, e.g., by interactions with high 
energy photons from the disk corona - e.g., Begelman, Blandford & Rees (1984) - or by 
proton-proton coUisions at parsec scale - Anyakoha, Okeke & Okoye (1988) -. 

Observations show that a number of jets are highly polarised. This is a further argument 
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in favour of jets made of e"^ pairs as there is no internal Faraday rotation or depolarisation 
(the rotation produced by the electrons is compensated by that of the positrons), while an 
ep-plasma gives rise to these effects. If jets were dominated by ep-plasmas, observational 
limits on the degree of Faraday rotation and depolarisation imply thermal electron densities 
< 10~^cm"^ - e.g., Walker, Benson & Unwin (1987) - or a minimum energy cut-off for 
the electrons fa 50MeV (Wardle 1977). This suggests a lack of thermal (cold) matter, at 
least, in the most polarised sources. However, as the upper bounds on the thermal electron 
density derived from the depolarisation argument depend on the magnetic field structure, 
they might be underestimated if there are field reversals within the VLBI beam. Finally, 
the detection of circular polarisation at parsec scales in several radio sources (Wardle et 
al.l998; Homan & Wardle 1999) suggests that, in general, extragalactic radio sources are 
mainly composed of e^-plasma. The argument is based on the fact that circular polarisation 
is produced by Faraday conversion, which requires that the energy distribution of the jet 
emitting particles extends to very low energies. This in turn indicates that e^-pairs are an 
important component of the jet plasma. 

Sikora & Madejski (2000) conclude that X-ray observations of blazars associated with 
OVV quasars impose strong constraints on the e^ pair content of jets of radio-loud quasars. 
According to these authors, pure electron-positron pair jets can be excluded because they 
may produce too much soft X-rays, while pure electron-proton jets may emit too few non- 
thermal X-ray radiation. Therefore, only jets may be viable where the number density of 
electron-positron pairs is much larger than the proton number density, but small enough for 
protons still being dynamically more important. Hirotani et al.(2000) find that parsec-scale 
jets cannot be dominated by baryons, if the electron density determined by the amount of 
synchrotron self- absorption necessary for an optically thick jet component, is of the same 
order as the one derived from the kinetic luminosity. 

For jets composed of a pair plasma additional questions arise. How can they maintain 
their stability, as for a given particle density they are lighter and hence less stable than jets 
made of an electron-proton plasma? Can they transport enough energy to fuel the observed 
kiloparsec scale radio lobes? Concerning the stability of electron-positron plasma fiows, Sol 
et al.(1989) demonstrate that the two-flow model is stable against excitation of electrostatic 
waves as long as the magnetic field is larger than a critical value Be- In their two-flow model, 
the relativistic beam is responsible for the VLBI jet and the observed supcrluminal motion, 
while the slower wind gives rise to the kiloparsec scale jet. Achatz & Schlickeiser (1993) 
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also consider the stability of e^-plasmas against the excitation of electromagnetic waves 
when the magnetic field exceeds Be- They find that the plasma may be stable over large 
distances provided plasma waves are damped thermally. The problem of the energy supply 
to large scales with a limited number of particles is turned around by Celotti (1998). She 
argues that the jets of low power radio sources mainly consist of e^-plasma, while those of 
powerful radio-loud quasars are made of ep-plasma. However, measurements of the circular 
polarisation of VLBA-jets of several FRII sources by Wardle et al.(1998) favour e^-jets also 
in powerful sources (at least, at parsec scales). 

Another important and not yet solved question concerns the impact of the composition 
on the morphology and dynamics of jets from sub-parsec to kiloparsec scales. Observations 
indicate that jets are relativistic at parsec scales - e.g., Laing (1996) -, then decelerate until 
they become sub-relativistic or mildly relativistic - e.g.. Bridle et al.(1994) - at kiloparsec 
scales with advance speeds of the terminal hotspots in the range 0.01c to 0.1c (Liu, Poo- 
ley & Riley 1992; Daly 1995), where c is the speed of light. At large scales the observed 
morphologies and deceleration of powerful radio sources is governed by interaction with the 
external medium. Several theoretical models have considered the gross features of kilopar- 
sec scale jets, i.e., their long term evolution. Begelman & Cioffi (1989) - BC89 hereafter - 
consider a simple model to describe the evolution of the cocoon. Self-similar expansion is 
suggested by Falle (1991). Komissarov & Falle (1998) explore the large-scale flow caused by 
classical and relativistic jets in a uniform external medium. They find that jets with finite 
initial opening angles are recoUimated by the high pressure in the cocoon and that the flow 
becomes approximately self-similar at large times. 

Numerical investigations have also addressed the kiloparsec scale regime. Two-dimensional 
Newtonian hydrodynamic simulations of axisymmetric light jets were performed by Cioffi 
& Blondin (1992) in order to understand the evolution of the cocoon, and were compared 
with the simple analytic theory of BC89. The grid resolution is quite high (15 zones per jet 
radius), but the simulations do only cover a relatively short period of the cocoon's evolution. 
The Newtonian simulations of Hooda, Mangalam & Wiita (1994) cover the evolution of ax- 
isymmetric extragalactic jets up to 10^ y. Their simulations include isothermal atmospheres 
with a density stratification given by a power-law, surrounded by an even hotter, but less 
dense intra-cluster medium where the jets accelerate and coUimate. Hooda & Wiita (1996) 
extended the results of Hooda et al.(1994) to three-dimensions, but their simulations cover 
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only a distance of 35 (initial) jet radii, and thus are too short to shed light on the long term 
evolution of real sources. 

Marti, Miiller & Ibanez (1998) - MMI98 hereafter - studied the long term evolution of 
powerful extragalactic jets on the basis of relativistic hydrodynamic simulations (up to an 
evolutionary time of 3 • 10^ y with a relatively low numerical resolution). The results are 
compared and interpreted with a simple generalisation of the model of BC89. They find 
an evolution divided into two epochs. After a transient initial stage, the jet's evolution is 
dominated by a strong deceleration. The jet advance speed becomes as small as 0.05c due 
to the degradation of the beam flow by means of internal shocks and the broadening of the 
beam cross section near the hotspot. 

The present work extends the investigations of MMI98 to much later evolutionary times 
and also addresses the question whether the content of thermal matter in powerful kiloparsec 
jets does influence their morphology, dynamics and emission properties. Like MMI98 we 
assume that the dynamics of the jets is dominated by the thermal plasma - e.g.,Sikora & 
Madejski (2000) - and, thus, a hydrodynamic approach is appropriate. For this purpose we 
have performed long-term simulations of axisymmetric, relativistic jets of extremely different 
composition, i.e., jets made of a pure leptonic or baryonic plasma, propagating into a uniform 
environment. 

2 MODELS 

2.1 Equation of state 

In almost all previous jet simulations an ideal gas equation of state (EoS) with constant 
adiabatic index 7 has been used. This is a good approximation in both the nonrelativistic 
(7 = 5/3) and the ultra-relativistic limit (7 = 4/3). However, when there exist very large 
temperature gradients in the flow, it is more accurate to use a EoS including a temperature 
dependent 7. For the simulations that we present below temperatures range from about 
10^ K in the ambient medium to 10^^ K in the hotspots. Thus, both relativistic and non- 
relativistic particles will participate on an equal footing. Additionally, extragalactic jets are 
likely composed of a mixture of particles of different masses, i.e., the value of 7 depends on 
the composition as protons become relativistic at higher temperatures than electrons. 

An equation of state that describes a mixture of ideal, relativistic Boltzmann gases has 
been derived by Synge (1957) (see also Komissarov & Falle 1998, appendix A). The Synge 
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EoS used in our simulations includes protons, electrons and positrons. The composition of 
the plasma only changes due to fluid mixing as the production or annihilation of electron- 
positron pairs can be neglected due to the low gas density in the kiloparsec scale (see, 
e.g., Ghisellini et al.l992). Assuming plasma neutrality, only one parameter is needed to fix 
the composition, e.g. , the mass fraction of the leptons Xi — (p^- + Pe+ ) / p where Pe- and Pe+ 
are the local rest-mass densities of electrons and positrons, respectively. Using the Synge 
EoS instead of a constant-7 EoS requires about 50% more computation time, because the 
iterative (Newton- Raphson) computation of T(£, p, Xi) {e being the specific internal energy) 
involves Bessel functions. 

2.2 The Pcirameter Space 

The following notation holds throughout the paper. Subscripts h and m refer to the beam 
and the external medium, respectively. The speed of light is set to c = 1. 

Assuming a uniform and static external medium a relativistic jet is fully described speci- 
fying, at a given inlet of radius i?^, the density contrast 77 between the beam and the ambient 
medium, the pressure ratio K — Pb/Pm between the beam and the ambient gas, the beam 
speed f b or equivalently the beam Lorentz factor Wh =1/ ^^(1 — f^), the beam Mach number 
Mh = Vb/csb (where Cgb is the beam sound speed), and parameters that depend on the EoS. 
For an ideal gas EoS only the adiabatic index 7 needs to be specified. In this case, a relativis- 
tic jet is freely scalable in size and density, but not in velocity due to the scale introduced by 
the speed of hght. Using the Synge EoS (§2.1) 7 is replaced by two parameters describing the 
chemical composition in the jet and in the external medium, e.g., the lepton mass fractions 
Xib and Xim, respectively. The Synge EoS also introduces an extra mass scale (due to the 
electron and proton masses), and hence an additional (mass) parameter, e.g., mo = PmRfj- 
With this set of parameters {f], K, Mb,Wb,mo, Xib, Xi^} a jet model is scalable under the 
transformations t ^ at,x ^ ax,p — > p/a^ (a = const). In the following this scale freedom 
will not be used, but pm and Rb will be fixed instead of mo. 

Current observations provide only constraints for the parameter space. The kinetic lu- 
minosity Lkin of a typical PR II jet is lO^^erg/s (Rawlings & Sunders 1991; Daly 1995). The 
initial jet propagation speed can be restricted by observations of CSOs (Compact Source 
Object) which have linear sizes of less than 0.5 kpc and which most likely represent a very 
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early evolutionary stage of a typical powerful radio source. Observed values of the hotspot 
propagation velocity are a; 0.2c (Owsianik & Conwayl998; Taylor et al.2000). 

The density, temperature and composition of the external medium are pm ~ 10~^mp/cm^ = 
1.67- 10~^'''g/cm^, lO^K and Xi^ — me/{me + mp) ^ 1/1837, respectively (e.g., Ferrari 

1998). Values of Ri, ~ O.Skpc can be inferred from observations of kiloparsec scale jets. But 
Rb is also constraint through other model parameters (see below). There are also observa- 
tional constraints on the likely value of the beam Lorentz factor at parsec scale - Wi, ~ 10, 
e.g., Ghisellini et al.(1993) -. Furthermore, measurements of the degree of polarisation as 
a function of frequency may be used to set an upper limit on the mean thermal electron 
density number which is Ue ~ 10"^ cm~^ (Perley, WiUis & Scott 1979; Burch 1979; Ghisellini 
et al.l992). 

We have fixed the values of 1/^^, Xi^, p^, and the initial propagation speed in all 
our models. The latter is set equal to 

Id a/% /^v 

Vj = —^—rVb , 1) 

with ?7r = PbhbW^ / Pmhm- This is the propagation velocity derived by Marti et al.(1997) for 
a pressure matched jet that moves in one dimension only (i.e., without sideways expansion). 
If the jet is not pressure matched, the propagation velocity can still be computed using 
the previous formula as long as the sound speed in the external medium fulfils Cgm << 1. 
The density contrast rj and the composition of the beam (Xn,) are not directly accessible to 
observations, i.e., these quantities are treated as free parameters. The beam Lorentz factor 
Wb cannot be chosen arbitrarily, because there are forbidden areas in the parameter space 
where the pressure ratio K and, most importantly, the beam temperature have unrealistic 
values (see below). We used values of Wb in the range 6.6 — 8. 

Lkin is obtained integrating the energy flux (see, Marti et al.l997, Eq. 20) over the beam 
cross section 

Lkin = {hbWb - l)pbWb7rRlvb . (2) 

This differs from the definition that other authors have used in the past - see e.g., Ghisellini 
(1998) - by the factor h^Wb — 1 which accounts for the fact that we do not include the 
rest-mass energy in the energy density, as it cannot be extracted from the beam particles. 

Having fixed Lkm, vj'^, Xib and the external medium, the values of rj, Th and M5 are 
uniquely determined by Wf, and Rf, through relations (1) and (2). To clarify this point, let 
us define function 
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C,:=,M^,T,) = ^^^^^-^ (3) 
and 

r?(Mr?, - 1) = ■ (4) 

Then = C^W^ — C^- For there exists no closed form, instead we have the relation: 

K(Cm-C.,T,) = ^^^^_^ (5) 

From Tfe and the EoS one can compute Csh and thus the Mach number Mb. Figure 1 shows 
the dependence of and r] on Wy, and R\y in the case of leptonic jets. For large values of 
W}) there is only a small range of values of around 0.35 kpc where solutions exist. In fact, 
the asymptotic value of the beam radius 



/5^.L_^ (6) 

is obtained when W;, — > oo and TJ, = 0. 

Computing Ri, from (6) using the independently measured values of Lkim vY and 
from observations yields T?^ ^ 0.35kpc, which is in agreement with standard values for kpc 
scale jets (e.g., Ferrari 1998). Fixing the values of and Lkim and for a constant large Wi, 
{i.e.^vi, ~ c), it turns out that 

k :— rjhb — const (7) 
and 

{kWb -r])Rl^ const . (8) 

The latter equation implies that for a given value of Wb, rj must increase with increasing 
it!;,. As k — const this means hb and thus Tb must decrease. Eventually, if Rb increases, Tb 
becomes negative. This explains the physically forbidden area to the right of i?^ in Fig. 1. 
A similar argument holds for the forbidden area to the left: decreasing Rb eventually leads 
to non-physical solutions with 77 < 0. Figure 1 also shows that for a given Wb there exists a 
maximum allowed r] and for a given r] there exists a maximum Wb- This maximum Wb grows 
increasing the value of v^. Let us also point out that from Eq. (4), increasing values of Lkin 
lead to larger values of Rb (notice the quadratic dependence). Hence, although there exits 
no solution with the seemingly reasonable values Wb — 10, rj — 10~^ and Rb ~ 0.35kpc (for 
the chosen parameters of the external medium, v^, and Lkin), one can obtain models with 
these characteristics by increasing and Lkjn up to 0.24 and 1.6 • lO^^erg/s, respectively. 
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Slightly smaller values of Wj, have been chosen to keep our parameters the closest possible 
to the observed values. 

Most of the previous jet simulations assumed that the jets are pressure matched {K — 1), 
because (i) there are no direct measurements of K, and (ii) this reduces the parameter space. 
However, fixing K implies that the properties of the external medium depend on the choice 
of the jet parameters. Therefore, one has to adjust the temperature of the external medium 
of every jet model such that K = 1. As wc wanted to simulate typical FRII-jets in a typical 
environment, the properties of the external medium are fixed and the assumption K = 1 is 
abandoned. For K > 1 jets can still remain stable, because they are pressure confined by the 
cocoon, where the pressure is much larger than in the external medium. In the simulated 
models K is in the range 1 — 200. Under-pressured jets {K < 1) cause considerable numerical 
problems when they are evolved beyond 3-10^ y, because the amount of ambient gas entrained 
into the beam becomes large. As the beams of our models are under-dense with respect to 
the external medium, dense blobs of matter from the ambient medium begin to pile up in 
front of the nozzle blocking the inflow. As we assume axisymmetry the material piles up 
around the jet axis. 

2.3 Numerical models 

The simulations are performed with the relativistic hydrodynamic, high-resolution shock- 
capturing code of MMI98. However, instead of a constant gamma-\sm EoS we have used 
the EoS described in Sect. 2.1. The code is well suited to solve the equations of relativistic 
hydrodynamics in cylindrical coordinates, i.e., for problems with axial symmetry. It has been 
extensively tested and apphed previously (e.g., , Marti et al.l995, 1997, MMI98). Covering a 
time span of up to 6.6 • 10^ y the set of models represent the longest simulations of relativistic 
jets performed so far. 

We have considered three models (Table 1) including two leptonic [Xn, — 1) ones and 
one baryonic model {X^, — 10~^; BC). One of the leptonic models has a factor of 100 
lower density 77 = 10~^ (LH) than the other two models with 77 = 10~^. The corresponding 
number densities of thermal electrons in the beam are ~ 2 • 10~^ and ~ 2 • 10~^, 
respectively. Note that both values are below the upper limit estimated from measurements 
of the degree of polarisation (see Sect. 2.2). The low 77 leptonic model has been run in order 
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Table 1. Parameters of the three simulated jet models. 



Model: 


BC 


LC 


LH 


Lkin [erg/sec] 




10*^ 




v}" [c] 




0.2 




Xih 


10-3 


1 


1 




10-3 


10-3 


10-5 


et [c2] 


5.69 • 10-3 


0.298 


119 


niK] 


lOio 


109 


2.37 • 10" 


Mb 


16.4 


2.38 


1.71 


Wb 


7.95 


6.62 


6.62 


K 


1.4 


91 


217 


lb 


1.42 


1.50 


1.33337 


Rb [kpc] 


0.366 


0.361 


0.342 



to check how extremely low values of rj affect the long term evolution of a relativistic jet 
(see, e.g., Birkinshaw 1991; Ferrari 1998). 

All models have the same power and initial propagation velocity, and they all propagate 
into the same homogeneous external medium with pm = 10~^mp/cm^, — lO^K, and 
= fT^e/ {f^e + f^p) ~ 1/1837. As thc jct powcr is fixed, the lower the density the higher is 
the internal energy, i.e., the light jet is also the hottest {e ^ 1). Note also that although the 
values of the injection temperatures (Table 1) seem to be rather large, they are unanimously 
determined once the fiducial conditions in the external medium, Lkim and Xu, have been 
fixed. 

We have not considered a light baryonic model (the counterpart of the light leptonic 
model) for three reasons. First, the temperatures within the beam and in some parts of 
the cocoon would reach 10^^ K. Thus, the dynamics of the model would be very similar to 
that of the hot leptonic model, i.e., radiation dominated. Second, at temperatures of 10^^ K 
(i.e., Lorentz factors of ^ 10^) leptons of the thermal plasma will contribute significantly 
to the synchrotron radiation power. Assuming an equipartition magnetic field (~ 60/xGauss 
for our simulations; see Sect. 5.4), the synchrotron cooling time is ^cooi ~ 10^ — 10^ y. This 
time is much smaller than the typical lifetime of the source, i.e., one has to consider syn- 
chrotron cooling in the simulations. Third, the light baryonic model is too hot. Its hot 
thermal plasma would have a detectable emissivity which is not observed in real sources, so 
far (see, e.g.,Celotti et al.l998). 

In order to track the evolution of the beam material the beam mass fraction X — Pb/ pis 
evolved by our code by means of an additional continuity equation. We also have to evolve 
the lepton fraction Xi because, as explained in Sect. 2.1, the composition may change due 
to mixing of species within the flow. 
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Table 2. Length and time units for each model. The third and fourth row give the total duration of simulation tmax in code 
and physical units, respectively. In the last row the time tcomp = 6.3 • 10® y is given in code units for each model. 



Jet 


BC 


LC 


LH 


Rb [kpc] 


0.366 


0.361 


0.342 


Rb/c [y] 


1192 


1176 


1114 


tmax [Rb/c] 


5300 


5400 


5950 


tmax [y] 


6.3 • 10*^ 


6.4 • 10*^ 


6.6 ■ 10® 


tcomp [^b/c] 


5300 


5350 


5650 



Our computational domain spans a region of size 200i?6x500i?ft (approximately 70 kpc x 175 
in cylindrical coordinates (r, z) with an uniform grid whose resolution is 6 cells per beam 
radius. This resolution is a compromise between a reasonable computing time per model 
(about 185 hours on a NEC SX-5 vector computer, running with a sustained performance 
of about 2.3 Gigaflops) and the maximum evolutionary time to be reached in a simulation 
(see Sect. 5.5). Axisymmetry is assumed along the r — boundary, while the downwind 
boundaries at r = 200i?b and z — 500Rb are outflow boundaries. The boundary at z = is 
reflecting. 

3 RESULTS 

3.1 Morphology and dynamics 

Due to the selected model parameters the beam radius and, hence, the time units are slightly 
different for each model. They are given in the first and second row of Table 2. The maximum 
evolutionary time reached in each model is listed in rows three (in code units) and four (in 
years), respectively. The last row gives the evolutionary time (in code units) at which all 
models are compared. This time corresponds to the final evolutionary time of the shortest 
simulation (model LH) and is equal to 6.3 • 10® y. 

Colour coded snapshots of the distributions of density, temperature and Lorentz factor, 
nearly a the end of the simulations, are displayed in Figs. 2-4. Additional contour lines mark 
the boundaries of the cocoon (see Sect. 3.3). At first glance, the morphology of all models is 
very similar, particularly concerning the overall cavity and beam shapes. 

Figure 5 displays profiles of several variables along the symmetry axis. One notices that 
the differences among models in quantities which are most important for the dynamics {p, 
e and P) are small, at least far from the head of the jet. The main differences occur where 
the profiles intersect the biconical shocks and rarefactions within the beam which are quite 
different (in spatial distribution and strength) from model to model. 
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The differences in the adiabatic index in the cavity (from model to model) are the result 
of the different initial beam temperatures and compositions. The leptonic models reach lower 
values of the adiabatic index in the cavity than the baryonic one (i.e., they have 7 closer to 
4/3) because electrons and positrons become relativistic at lower temperatures than protons. 
Although model LC has an average temperature within the cavity (see Fig. 3) more than 
ten times smaller than in model BC, as a result of the mixing with the ambient medium, 
the effective adiabatic index in the cavity is smaller for model LC than for model BC. The 
reason being that model LC supplies with a thousand times more relativistic particles the 
cavity than model BC. The beam of model LH is the hottest and, thus, it has the lowest 
adiabatic index. It is even lower than model LC because the leptons of model LH are more 
relativistic (i.e., the beam temperature is higher). 

At the jet head different phases of vortex shedding can be seen in Fig. 2. In model LC 
new vortices are forming, in model BC they have been just shed from the head, and in model 
LH the Mach disk has been replaced by a conical shock causing an acceleration of the head. 
However, these differences are present only temporarily. More important are the differences 
inside the beam and in the cocoon which we discuss next. 

3.2 Beam 

In the hot, light model LH the density in the beam is two orders of magnitude lower and the 
internal energy is, at least, two orders of magnitude larger than in the two colder, denser 
models BC and LC. However, model LH has a strongly pinched structure (at z — 39 Rb and 
z — 63Rb in Fig. 4) while the colder models possess a continuous beam where the beam mass 
fraction is larger than 0.95 throughout the beam. The beam pinching, even disruption, is 
due to considerable mass entrainment in model LH, inside which one can find regions which 
are almost at rest. Hence, the light, hot jet appears to be the most unstable one. 

The thermodynamic quantities change drastically at the first conical shock in model BC. 
The temperature rises from lO^^K to lO^^K, the pressure (initially model BC is pressure 
matched, i.e.,K fa 1) increases by two orders of magnitude and stays at values similar to 
those of the other models from there on {K — (9(100)). The internal energy and sound 
speed also strongly increase, while the relativistic Mach number drops from 130 to 20. All 
these changes are caused by the initially very low internal energy of model BC (see Tab. 1). 
The other two models initially have a much larger internal energy and the thermodynamic 
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quantities in the beam change less during the evolution. Actually, the properties of the 
recoUimation shocks are strongly varying from model to model. On average, recoUimations 
shocks have larger compression ratios and are more numerous in model LH than in the colder 
models (see Fig. 5a). 

The sound speed approaches its maximum value (c^ ~ 0.57) close to the injection nozzle 
in model LH. In the colder models Cg increases within the beam and the maximum value 
is reached at the hot spot. The deceleration (acceleration) at internal shocks (rarefactions) 
within the beam is much more violent in model LH (where Lorentz factors of up to 20 can 
be observed) than in the two colder models. 

At the beam/cocoon interface a thin, hot layer forms in the cold models, which is thin- 
ner than that found in 3D simulations (Aloy et al.l999, 2000). This layer results from the 
interaction of the beam with the external medium, and appears naturally in some models of 
jet formation (Sol et al.l989). Its existence has been proposed by different authors (Komis- 
sarov 1990; Laing 1996; Laing et al.l999) in order to account for a number of observational 
characteristics of FRI radio sources. However, the physical nature of the shear layer is still 
largely unknown, as a study of its properties requires much better resolution and full 3D 
simulations. The structure of the shear layer is remarkably different in model LH, where the 
hot layer forms out of the beam/cocoon interface due to the extremely high specific internal 
energy of this model. Nevertheless, the variations of pressure and Lorentz factor across the 
shear layer are too small in model LH to have an effect on the emission. This is in contrast to 
Aloy et al.(2000), where the shear layer was much more extended and its emission properties 
were distinguishable. 

Despite of the differences between models, all of them maintain relativistic beams up 
to distances of more than 80 kpc the average Lorentz factor being larger than 5. This has 

important theoretical consequences, in particular, for the most commonly accepted emission 
models of kiloparsec-scale X-ray jets as we will discuss in Sect 5. 

3.3 Cocoon 

The cocoon surrounding the beam is formed by beam matter that flows backward (relative 
to the beam) after being deflected at the hot spot. This beam matter partially mixes with 
the external medium, i.e., the cocoon does not contain pure beam material {X ^ 1). Hence, 
we define the cocoon as the region containing material with a beam particle fraction 0.1 < 
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X < 0.9. We have used the beam mass fraction as the criterion, because the density within 
the cavity is similar in all models (see Fig. 2), i.e., the lower the value of X, the lower is the 
number of emitting particles and thus the total emitted intensity. The thresholds are a bit 
arbitrary, but we have checked that variations of the lower threshold (up to a decade) do not 
change significantly the shape of the cocoon. For a very small lower threshold {X^nin — 0.001) 
the cocoon practically coincides with the complete cavity blown by the jet (see below). 
Increasing the value of the upper threshold would include too much of the beam itself as 
part of the cocoon. Other cocoon definitions have been given in, e.g., Cioffi & Blondin (1992). 

According to our definition, the shape of the cocoon is remarkably different among the 
considered models. In model LH the cocoon is restricted to a narrow layer around the beam, 
while it is much more extended and matches the relativistic backfiow region in the denser 
models (BC and LC). This is due to the lower beam density of model LH, which allows for a 
larger inertial confinement of the beam fiow by the external medium without causing much 
mass entrainment (which would destroy the jet). Another consequence of the lower density 
of model LH is a lower beam mass fraction within the cavity outside the cocoon (about two 
orders of magnitude smaller than in the other models). The average cocoon temperature of 
models BC and LH is around lO^^K the size of local fluctuations reaching up to two orders 
of magnitude. In model LC the average cocoon temperature is lower (10^°K) and there are 
almost no fluctuations (see Fig. 3). The temperature distribution of models BC and LH is 
very similar (although they have a different density and composition), whereas that of model 
LC (with the same density as model BC and the same composition as model LH) is very 
different. In Sect. 5 we will explain how the varying particle density causes this effect. 

A backfiow with mildly relativistic velocities occurs in thin shells starting from the head 
of the jet and limiting the cocoon in models BC and LH. The amount of beam material in 
this backfiow is smaller in model LC than in the other models (see, e.g.. Fig. 4). 

4 EVOLUTION 
4.1 Definitions 

For the discussion of the jet evolution we introduce the following quantities: (1) the jet 
length Ij being the z-coordinate of the contact discontinuity; (2) the cocoon length Icc being 
the difference between Ij and the smallest z-coordinate of the fiuid belonging to the cocoon 
(according to the definition given in Sect. 3.3); (3) the cavity length 1^ being the z-coordinate 
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of the bow shock on the jet axis; (4) the average cocoon radius r^; (5) the average cavity 
radius Vc, (6) the cocoon aspect ratio Aj = Icc/vj] (7) the cavity aspect ratio Ac = Ic/rc] 

(8) the hotspot pressure Phs being the average pressure in the 10 cells upstream of 2; = If, 

(9) the average cavity pressure Pc, (10) the jet velocity Vj being the velocity of the contact 
discontinuity between the jet and the external medium; and (11) the bow shock velocity Vc 
on the jet axis. 

The position of the contact discontinuity (Ij ) is always near the bow shock (/c, see Fig. 6a). 
But whereas lc{t) is increasing monotonically, lj{t) looks rather 'noisy' due to the vortex shed- 
ding mechanism, which makes it difficult to detect (numerically) the position of the contact 
discontinuity. The detection algorithm searches for the first cell on the jet axis (beginning 
at 2; = 0) where X < 0.5. When a vortex sheds off, the distance of this cell from the inlet 
is abruptly reduced. As Aj and Vj depend on Ij, they display an oscillatory behaviour, too 
(Figs. 6b and 6d, respectively). Other quantities like rj and Phs do also oscillate because 
of the highly dynamical processes that affect them (vortex shedding and mixing with the 
external medium), in particular when the jets enter their second evolutionary phase (see 
Sect. 4.2). Pfis and Vj have been smoothed (averaged over 25 and 50 cells, respectively) in 
Figs. 6a and 6c, but their almost constant values during the first 100 time units are not 
caused by this smoothing process. 

4.2 Evolutionary phases 

In all models two evolutionary epochs can be distinguished. During the first epoch {t < t^^ ^ 
1.2 • lO^y; see Fig. 6) multidimensional effects are not yet important and the jet propagates 
with a velocity close to the ID estimate (1). The Mach disk is only shghtly disturbed, the 
shedding of vortices is negligible, and both the hotspot pressure and the jet velocity are 
nearly constant (Figs. 6c, 6d). The aspect ratios of the cocoon and of the cavity (Fig. 6b) 
increase with time, while the average cavity pressure decreases (Fig. 6c). 

The second epoch starts when the first larger vortices are produced near the jet head. 
The beam cross section increases and the jet decelerates progressively. The deceleration is 
repeatedly interrupted by short acceleration phases caused by the temporary replacement 
of the Mach disk by a conical shock, which occurs e.g., in model BC at i ~ 400, 1200 and 
5100, respectively (Fig. 7). The velocity increase during the acceleration phases can exceed 
100%. 
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Table 3. Cavity lengths (in kpc) and aspect ratios at the end of the one-dimensional epoch (t = 1.1 • 10^ y; first two rows) and 
near the end of the simulations (t = 6.0 • 10® y; third and fourth rows) for all models. The last row shows the time averaged 
velocity during the time interval [1.1 • 10^,6 • 10®] (in units of c). 



Phase- 


.Id 


!!(' 


LC 


LH 


ID 




7.57 


7.10 


6.83 




Ac 


3.37 


3.09 


2.89 




Ic 


106.8 


100.6 


97.4 


2D 


Ac 


4.32 


4.03 


3.79 




V 


0.055 


0.052 


0.050 



The acceleration phases also give rise to a higher hotspot pressure. Such pressure can 
be estimated assuming that the pressure of the external medium is negligible. Then the 
ram pressure of the external medium Pmvf. is equal to the thermal pressure behind the bow 
shock, and hence about equal to the hotspot pressure as the pressure is continuous across the 
contact discontinuity (Fig. 8). There is a phase shift between P^g and Pm'^l, because pressure 
fluctuations are generated at the contact discontinuity (vortex shedding) which arrive at the 
bow shock later. 

5 DISCUSSION 

In spite of an extremely different composition and internal energy, the differences in mor- 
phology and dynamics of the jets are small (particularly, between models BC and LC, both 
having the same density contrast). This is unexpected, as the jets of our models have a fixed 
power which is carried by particles of different mass, i.e., the energies per particle are quite 
different. A leptonic jet has fa 10^ times more particles than a baryonic jet of the same mass 
density. Therefore, its kinetic energy per particle is smaller, i.e., its temperature is smaller, 
but its specific internal energy (erg/g) is larger. We also expected that differences in the 
adiabatic index of the models would manifest itself. 

It turns out, however, that fixing the power and initial speed of the jet as well as the 
properties of the ambient medium, the development of the jet seems to be quite well defined 
(at least for the jet parameters considered and for the evolution time covered by our sim- 
ulations). However, a more careful inspection of the results reveals that some remarkable 
differences exist between our models which we discuss in the following. 
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5.1 Comparison with the extended Begelman-Cioffi model 

A simple theory of the evolution of the cocoon of extragalactic jets was presented in BC89. 
The authors argue that the cocoon has not yet reached pressure equilibrium with the sur- 
rounding medium in many sources, and that high pressure confines the jet keeping it highly 
coUimatcd as seen in extragalactic sources. BC89 assume that both the bow shock velocity Vc 
and the power transported into the cavity are not time-dependent, and that the pressure 
of the external medium is negligible. Then the average cavity pressure is given by 

where Fc — vrr^ is the cross section of the cavity and 7c = Cp/c^ is the constant average 
adiabatic index in the cavity. The cavity pressure causes an expansion of the cavity with a 
velocity fc, i.e., 

Pc = Pmf-l . (10) 

This implies l/rccafc and hence 

Tc OC t^/^ Fc (Xt, Pc(X Ij/Vc OC t^/^ . (11) 

We have extended the model of BC89 (eBC, hereafter) replacing the assumption Vc — 
const by the more realistic one Vc oc t"' and further assume that fc (X t^. Then (9) together 
with (10) yields 

^ /3 = -l/2-a/4, (12) 



and thus 

rcoct^ Fc(xt 2^ PcOct 2^ Ij/vc (X t2^ ^ . (13) 

The only free parameter of the eBC model is the value of a. A deceleration (a < 0) leads to 

a faster radial expansion, a slower decrease of the cocoon pressure and to a slower increase 
of the aspect ratio. For a = —2/5 the cavity evolves self-similarly (i.e., with constant aspect 
ratio) . 

We have extracted the parameter a and the exponential dependence of Pc and Ij/vc from 
our simulation results by fitting the position of the bow shock as a function of time with 
a power law (Tab. 4). Using the fitted values of a we have also computed the exponents 
for P°^^ and (/j/Tc)°^^ according to the extended BC89 model. The fitting procedure is 
not very reliable for the ID epoch as it involves only a few data points, while velocity 
fluctuations cause problems for fitting the 2D epoch. Thus, differences between the models 
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Table 4. Comparison of the simulation results with the extended BC89 model. Columns one and two give the evolutionary 
epoch and the model, respectively, and column three shows the value of the free parameter a obtained from the simulations by 
fitting the position of the bow shock as a function of time with a power law. The fourth column provides the fitted exponential 
time dependence of the cavity pressure, while the fifth column gives corresponding exponent of the eBC model (= — 1 — a/2). 
Columns 6 and 7 show the exponent describing the time dependence of the aspect ratio obtained from the simulations and the 
eBC model(l/2 + 5a/4), respectively. 



Phase 


Jet 
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Pc 


pcBC 
^ c 


Ij/r-c 


(/,/r,)=BC 




BC 


-0.113 


-0.945 


-0.944 


0.449 


0.358 


ID 


LC 


-0.151 


-0.940 


-0.925 


0.404 


0.311 




LH 


-0.219 


-0.906 


-0.891 


0.358 


0.226 




BC 


-0.355 


-0.722 


-0.823 


0.058 


0.056 


2D 


LC 


-0.356 


-0.737 


-0.822 


0.071 


0.054 




LH 


-0.363 


-0.779 


-0.818 


0.059 


0.045 



are probably not significant. Note in this respect that although the value of a varies by a 
factor of two between the models, the variation of P^^'^ is at most 5% and that of {Ij/tcY^^ 
less than 25% as both quantities can be extracted more reliably. The values of P^^*^ and 
{Ij/vcY^'^ obtained from the simulations and from a via the eBC model agree reasonably 
well. The exponent for the aspect ratio is shghtly larger than that predicted by the eBC 
model (especially for the ID epoch), i.e., the radial expansion of the cavity is slower than 
expected from the deceleration of the jet head. This discrepancy arises from the assumption 
of the BC89 model that the whole power of the jet is used to increase the pressure behind 
the bow shock (which is responsible for the sideways expansion of the cavity). However, part 
of the jet power is used to fuel the relativistic backflow and its kinetic energy is not available 
to expand the cavity in radial direction, i.e., the aspect ratio (jet length over cavity radius) 
grows faster in the simulations than in the BC89 model. 

During the 2D epoch the cavity evolution is almost self-similar and the cavity aspect ratio 
remains nearly constant. This situation will only change when the cavity pressure becomes 
equal to the pressure of the external medium and the radial expansion ends. Extrapolating 
from Fig. 6 this will happen after ~ 10^ time units or 10^ y, which exceeds the estimated 
ages of observed jets (10^ y - 10^ y) by one order of magnitude. 

The evolution of the jets proceeds in two epochs independent of the composition of the 
jet. This is in agreement with previous results obtained by MMI98 and Komissarov & Falle 
(1998). The ID epoch of our models corresponds to the intermediate phase of Komissarov & 
Falle (1998), while the second epoch may be identified with their self-similar phase. They also 
find that the growth rate of the aspect ratio is smaller during the second epoch. However, 
their growth rates are slightly larger than ours in both epochs. This discrepancy is probably 
due to different jet parameters. Their relativistic models are 30 times denser than ours, have 



© 2001 RAS, MNRAS 000, 1-?? 



Does the plasma composition affect the long term evolution of relativistic jets? 19 
a little bit smaller Lorentz factor, are pressure matched, and possess a finite opening angle. 
The transition to the second epoch does not occur in those models of Komissarov & Falle 
(1998) which have an opening angle larger than 10°. 

5.2 Evolutionary differences 

The average jet velocities during the 2D epoch are listed in Tab. 3 together with the lengths 
and the cavity aspect ratios at the end of the ID epoch and at the end of the simulation. 
At first glance there seems to be a clear trend in the ID epoch: with increasing internal 
energy and decreasing Mach number the jets become slower and wider (the mass flux into 
the cocoon increases). As already found by MMI98, the jet velocities are noticeably larger 
than predicted by the ID estimates (wj'^ = 0.2c) at the beginning of the simulation (see 
Figs. 6d and 7). This discrepancy is caused by imposing a reflecting boundary condition at 
z = (see Sect. 2.3) and by the different pressure ratios K. After some ten time units the jet 
head has propagated far enough from the grid boundary to be no longer affected by it. The 
average cavity pressure decreases and the pressure contrast between the jet and the cavity 
becomes closer to one in all the models (Fig. 6). After the jets has decelerated at the end of 
the ID epoch their velocities and the ID estimates agree very well (Fig. 7). 

During the 2D epoch the jet of model BC propagates considerably faster than those of 
models LC and LH, while the jet of model LC is only slightly faster than that of model LH. 
This result can be easily understood. The relativistic Mach number of model BC is larger, 
i.e., the angle between the conical shocks inside the beam and the beam axis is smaller. Hence, 
these shocks are less efficient in decelerating the flow. This also affects the propagation speed 
of the head, because the terminal Mach disk is temporarily replaced by conical shocks an 
effect pointed out already by Marti et al.(1997). The noisy bow shock velocity (Fig. 7) is 
also a result of this process. However, as the heads of all jet models are hot and the Mach 
number is low in this region, this effect does not lead to very strong differences between the 
computed models. 

The differences in the cavity aspect ratios correlate with the average propagation speeds. 
The jet of model BC is the fastest and the longest one (Tab. 3). The cavity radius is similar 
for the cold models and bigger for model LH. The models expand radially almost at the same 
rate, because the average cavity pressure is essentially equal for all three models although a 
signiflcant part of the kinetic energy which is contained in the relativistic backflow cannot 
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be used to inflate the cavity. In fact, we included model LH in our investigation in order 
to check whether an increased fraction of internal (non-directed) energy would lead to less 
kinetic (directed) energy in the cocoon. 

5.3 The influence of different compositions 

From a numerical point of view the most important effect of a varying composition is the non- 
constant adiabatic index requiring an extension of the original MMI98 code (see Sect. 2.3). 
The composition dependent 7 calculated from the Synge EoS (see Sect. 2.1) however, does 
not lead to large differences in the evolution and morphology of the jets. The average 7 in 
the cocoon is around 1.4 and 7 = 5/3 in the external medium for all three models. Due to 
the initial conditions, the adiabatic index within the beam is different for every model, but 
these differences have no observational relevance as 7 is not an observable quantity. 

We expected that the largest internal energy content in the leptonic models would cause 
a larger expansion of the beam than in the baryonic model, leading to a faster deceleration 
of the head and differences in the aspect ratio of the cavity (more spherical in the leptonic 
models). However, this effect is compensated by the large overpressure of the cavity that 
efficiently confines the jet laterally. 

The variation of the particle density with the composition gives rise to huge differences 
in the temperature. At the same density, an e^-gas contains three orders of magnitude more 
particles than an ionised hydrogen gas. The particle density is proportional to rj ■ Xn, and, 
therefore, it is much higher in model LC than in the other two models. Hence, the same 
amount of energy is distributed over many more beam particles in model LC than in models 
BC and LH, leading to an about three orders of magnitude lower beam temperature (Fig. 3). 
However, this does neither explain the fiat temperature profile of model LC, nor the identical 
average cavity temperatures of models LH and BC (whose hotspot temperatures differ by 
a factor of ten). To explain the temperature differences in the cocoon, particles from the 
external medium must be considered too. There are A^cm — 7- 10^^ particles from the external 
medium and Neb — 3-10^^ particles from the beam in the cavity of model LC ait — 6.0- 10^- 
The dominant fraction of the internal energy of the cavity is also carried by beam particles. 
Hence, the mixing of beam and external medium has little impact on the internal energy 
per particle, i.e., the temperature remains at the hotspot value. 

Concerning models LH and BC, particles from the external medium dominate the particle 



© 2001 RAS, MNRAS 000, 1-?? 



Does the plasma composition affect the long term evolution of relativistic jets? 21 
density in the region of the cavity excluding the dense shell behind the bow shock {cavity 
interior), i.e., Ncib/Nc^m < 1, but the internal energy is still dominated by beam particles. 
This leads to large changes and fluctuations of the temperature during the mixing of particles 
from the beam and the external medium (Fig. 3). The temperatures near the jet head and 
near the inlet are very different in these two models. Although the hotspot temperatures 
differ, the average cocoon temperatures are similar because the same amount of energy is 
distributed over a similar number of particles. In addition, since the particle density is much 
lower than in model LC, the average cocoon temperature is much higher. 

In general, when Nc^b/Nc^m > 1 the cavity interior has a flat temperature profile and 
we talk of an isothermal cavity. When Nc^h/Nc^m < 1 the average temperature in the cavity 
interior is heterogeneous and independent of the beam composition and density. In this case 
we talk of a non-isothermal cscvity. Another effect is caused by the different time dependencies 
of Nc^b and Nc^m- The number of particles from the external medium in the cavity interior 
is proportional to its volume, and the cavity interior has a roughly constant density during 
the whole evolution (e.g., , it changes by less than 20% during 6.0 • lO^y in case of model 
LH). Hence, from the eBC model one can derive A^c^m oc i2(i+/3) . ^i+a ^ ^2+^/2 while the 
number of beam particles increases in the cavity interior linearly with time {N(.^b oc t) , the 
ratio Nc^b/Nc^m oc t~^~°'/^ decreases, i.e., a transition from an isothermal to a non-isothermal 
cavity may occur during the jet evolution. 

Figure 9 shows the evolution of Ncib/N^m for the three models. Model LC has an isother- 
mal cavity throughout the simulation, while the cavities of the other two models are always 
non-isothermal. Although, we have not found the above mentioned transition. Fig. 9 suggests 
that one should be able to find initial jet parameters (between those of models LH and LC) 
where such a transition will occur during the simulation. 

In an isothermal cavity the external particles can be neglected, i.e., Nc ~ Nc^b oc t. As the 
cavity volume is proportional to the number density within the cavity Uc oc 

However, as Pc oc this implies that the temperature Tc oc Pc/uc does not depend 

on time. This explains why the cocoon temperature is constant in 2;-direction and also in 
r-direction (Fig. 3). In a non- isothermal cavity particles from the external medium dominate, 
i.e., Uc — const and the average cocoon temperature decreases according to Tc oc t"^""/^. 
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5.4 Comparison with observations 

Comparing the hydrodynamic properties at the end of our simulations with those deduced 
from observations of CygA, we find that the hot spot pressure {Phs ~ 0.001p,„c^ ^ 1.5 • 
10~^dyn cm~^) agrees within a factor of 2 with the value reported by Carilli et al.(1996). 
Comparing the hotspot advance speeds of CSO and FR II sources implies a deceleration of 
the hotspot propagation at early epochs. This deceleration indeed occurs in our simulations. 
The propagation velocity after the deceleration {vj ~ 0.05c) is in good agreement with 
the estimates of Daly (1995). The equipartition magnetic field of ~ 60 //G in the cavity 
and ~ 600 //G in the hotspot are larger than those reported for Cyg A (Carilli et al.l996), 
although they agree with the values obtained for the radio lobes and hotspots of younger 
powerful radio sources (Ferrari 1998). This agreement is easily explained considering that 
the evolutionary time covered by the simulations is relatively short. 

Up to now it is impossible to determine hydrodynamic properties of jets directly from the 
observed radio emission. Thus, in order to compare the simulated jets with observed ones 
we proceed in 'opposite' direction and compute the non-thermal synchrotron radio emission 
of our hydrodynamic models. We follow the same approach as in Aloy et al. (2000) assuming 
that the energy density weighted with the beam particle fraction is distributed among the 
emitting non-thermal particles according to a power law N{E) oc E^'^ (without cut-offs). 
The magnetic field is considered to have equipartition strength and a ad hoc structure: 
aligned with the flow velocity and with a negligible random component (because we are 
not concerned in polarisation properties). The simulated radio maps obtained from the 
hydrodynamic models (Fig. 10) include only those computational zones where the density is 
smaller than O.lpm to exclude emission from the external medium. We have chosen a viewing 
angle of 45° to compute the synthetic radio maps in order to avoid excessive Doppler boosting 
of the beam which would out-shine completely the diffuse emission from the cocoon. 

Two main features known from observations are present in all models: radio lobes with 
hotspots and a one-sided, knotty jet (particularly evident in model LH). The knots are 
associated with internal shocks in the beam. Emission is dominated by the hotspots, where 
the pressure is maximum. Figure 10 shows that the models have different jet/cocoon emission 
ratios, as expected from their different cocoon morphology (Sect. 3.3). While models BC 
and LC show some emission from the cocoon, resembling the observed lobes, the emission 
is dominated by the beam in model LH. 



© 2001 RAS, MNRAS 000, 1-?? 



Does the plasma composition affect the long term evolution of relativistic jets? 23 
In actual sources the lobes of the jet and counter-jet are often observed to be equally 
prominent. Our simulations show larger lobes (i.e., more cocoon emission) for the counter- 
jet (particularly for models BC and LC). This can be explained by the relativistic backflow 
present in the simulated models, because Doppler boosting enhances the emission from 
the backflow regions in the counter-jet and dims that of the approaching jet lobe. The 
relativistic backflow also limits the radial expansion of the lobes and prevents the beam from 
inflating large lobes with gas moving at sub-relativistic speeds. Therefore, the similarity of 
the lobes in actual radio sources suggests that there are no relativistic backflow regions. 
Their presence in the simulated models is most likely due to the axial symmetry imposed 
in our 2D simulations. Three dimensional simulations (Aloy et al.l999) show much smaller 
backflow velocities, typically ~ 0.25c, and hence do not lead to substantial Doppler beaming 
in the counter-jet cocoon. 

Very recently, the new generation of space-based X-ray observatories has allowed for 
the detection of X-ray emission from kiloparsec scale jets (Schwarz et al.2000; Chartas et 
al.2000; Sambruna et al.2001). The observed X-spectra are not satisfactorily explained by 
standard radiation mechanisms, i.e., synchrotron self-Compton (Tavecchio, Maraschi & Ghis- 
ellini 1998; Schwarz et al.2000; Tavecchio et al.2000) and synchrotron radiation. The latter 
can arise from the same population of particles that produce the radio emission (Schwarz 
et al.2000; Marshall et al.2001), or from a second, much more energetic population of elec- 
trons co-spatial with the one responsible for the radio emission (Roser et al.2000; Sambruna 
et al.2001). Alternatively, inverse Compton scattering of beamed photons of the cosmic mi- 
crowave background (CMB) has also been proposed (Tavecchio et al.2000; Celotti et al.2001). 
For this mechanism to work, it is necessary that the beam of the jet is, at least, mildly rel- 
ativistic, because the amphfication of the CMB radiation increases with the Doppler factor, 
5 = [r(l — (3cos9)]^^ (Tavecchio et al.2000). Our models show relativistic beam velocities 
(r ~ 5) out to 80 kpc (see Sect. 3.2). Whether this is a result of the imposed axisymmetry of 
our models or a physical feature may only be disentangled by performing three-dimensional 
simulations (see Sect 5.5). 

Thermal bremsstrahlung from the gas confining the radio-optical jet was very soon dis- 
regarded as the source of the observed X-ray radiation from the jets of powerful sources 
(Sambruna et al.2001; Schwarz et al.2000), although the data can be fitted by a thermal 
bremsstrahlung spectrum (Chartas et al.2000). In order to explain the observations invoking 
thermal bremsstrahlung, a large electron density {ug ~ 2 cm~^) is required. This is, however, 
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inconsistent with the observed (relatively low) rotation measure {ue < 3.7 x 10~^/{B L); B 
and L are the magnetic field strength and the path length, respectively; Schwarz et al.2000). 
The estimate for the electron number density, is based on the assumption that (i) there are 
no magnetic field reversals within the VLA beam, and (ii) the confining thermal material 
is not relativistic and relatively cold. But, our simulations show the presence of a very hot 
relativistic plasma (at least in model BC), which is not accurately described by the classical 
formulas. 

Given that the densities are very similar in all our models and that there arc huge tem- 
perature differences inside the cavity (Fig. 3), we have analysed whether relativistic thermal 
bremsstrahlung can be used to distinguish between the cavities infiated by jets of different 
composition. We follow the work of Nozawa, Itoh & Kohyama (1998), which includes the 
appropriate Elwert (1939) factor to compute the relativistic thermal bremsstrahlung cross 
section. Their method is accurate if one can neglect the thermal motion of protons, i.e., at 
temperatures below 10^^ K. For larger temperatures the relativistic thermal bremsstrahlung 
cross section is still accurate within an order of magnitude (Itoh, private communica- 
tion). Above 10^ K electron-electron, electron-positron, positron- ion, and positron-positron 
bremsstrahlung processes are important (see, e.g., Novikov & Thorne 1973). These processes 
should have been included when computing the total bremsstrahlung emissivity of our mod- 
els, where temperatures are as high as 10^^ — 10^^ K in the shear layer confining the jet, in 
the hotspot, and in the cocoon of model BC. However, as the emissivity of these processes 
is comparable with that produced by the interaction of electrons with protons (Dermer 
1986; Svenson 1982), the emissivity of the latter process provides a lower bound of the total 
bremsstrahlung emissivity, which is accurate to an order of magnitude up to temperatures 
of ~ 10^3 K 

The bremsstrahlung power p^''^™-^ per unit of frequency (z/) in the comoving frame is 
proportional to neni/\/T {rii is the number density of ions in the medium). Considering 
a viewing angle of 90°, i.e., beaming effects are unimportant, P^^'^^ is dominated in all 
models at relatively low frequencies (/ii/^lOOkeV) by the very dense shell at the cavity 
boundary behind the bow shock containing shocked external medium. Although the non- 
shocked external medium is cooler than the shell, P^^'^^ is three orders of magnitude smaller 
than in the shell, because the shocked material is much denser than the ambient medium. 
The bremsstrahlung emission from the cavity's interior is very weak in models BC and 
LH. This is different for model LC, where the temperatures in the cocoon are lower and 
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the electron densities are larger than in the other models. The resulting bremsstrahlung 
emissivity of the cavity's interior is almost as strong as that of the external medium. 

The dominance of the X-ray emission from the shell is also evident from the surface 
brightness, i.e., from the integral of p^^'^^ computed along the hue of sight (Fig. 11). The 
total emission X-ray maps of all models look very similar, as the shell outshines the cavity 
at a frequency of lOkeV (Fig. 11 top panel). Note that our analysis does not take into ac- 
count the emission of intra-clustcr gas, which would be dominant for a radio source located 
in a galaxy cluster. The dominance of the shell is in agreement with the work of Heinz, 
Reynolds & Begelman (1998). They also find that the shell's dominance decreases as the 
source evolves because the bow shock decreases in strength when the pressure inside the 
cavity decreases due to the expansion. McNamara et al.(2000) find a reduced X-ray sur- 
face brightness when observing the radio lobes of Hydra A, an effect which is evident in 
our models. However, McNamara et al.(2000) also point out that there is no evidence for 
shock-heated gas surrounding the radio lobes, and hence suggest that the cavity expands 
subsonically. This does not contradict our results, as we consider FRII sources (Hydra A is a 
FR I source) and because our models represent a still relatively young evolutionary stage of a 
powerful radio source. FRI sources reach the transonic regime relatively soon and, therefore, 
one does not expect to find a strong shock surrounding the radio lobes. Instead the radio 
lobes fiare into the external medium. This may explain the lack of X-ray emission from the 
region surrounding the cavity. On the other hand, it is expected that a further evolution of 
the source will lead to a decrease of the cavity pressure, and consequently to a weakening of 
the cavity bow shock (which in turn reduces the X-ray emission of the shell) . 

We have also computed the surface brightness per unit of frequency at lOMeV (Fig. 11 
bottom panel). We choose this high frequency, which is beyond the observational range of 
CHANDRA (~ 0.1 — lOkeV), because only when considering hard X-rays one begins to 
see the emitting fluid inside the cavity and the emission of the model differs most. Increas- 
ing the frequency from lOkeV to lOOkeV^i/^1 MeV the interior of the cavity shows up, 
although the bow shock surrounding the head of the jet is still dominant. At even higher 
frequencies (^1 MeV), the hotspot and the high temperature cocoon are visible. Finally, at 
hv — 10 MeV (bottom panel in Fig. 11), only models BC and LH show some emission while 
the exponential cut-off of the bremsstrahlung spectrum reduces the emissivity of model LC, 
as the temperature of the cavity interior of model LC, kT ~ 500 keV, is much smaller than 
the considered frequency. This might provide a way to distinguish observationally - once fu- 
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ture detectors have the appropriate dynamical range and sensitivity - leptonic from baryonic 
jets (with the same kinetic power and density contrast). If one observes a transition from a 
shell dominated source to a source dominated by the interior of the cavity when increasing 
the observation frequency from hu ^ 1 MeV to hu 10 MeV, the jet feeding the radio lobe 
is most probably rich in baryons. Conversely, if this transition is not detected and one finds 
instead a significant suppression of the emission, the jet is mostly made of leptons. 

Typical luminosities of very powerful X-ray, large-scale jets are of the order 10^^ erg sec^^, 
e.g.,Chartas et al.(2000) report a X-ray luminosity in the 2 — lOkeV band of Lx — 6.3 • 
lO^^ergsec"^ for PKS 0637-752. The X-ray luminosities produced by our jet models due to 
thermal bremsstrahlung are far below the observed values. This rules out thermal bremsstrahlung 
as the dominant mechanism for the X-ray emission of large-scale jets. However, several fac- 
tors may enhance the total thermal luminosity. First, the size of our jets is still relatively 
small as compared with fully evolved FR II sources (the larger the jet the larger is the lu- 
minosity). Second, the kinetic power of our models, although being large, is still below the 
one inferred e.g., for PKS 0637-752 (10^® erg see" ^; Tavecchio et al.2000). Our simulations 
show that the jet itself radiates only a small fraction of the total power per unit of frequency 
emitted by the whole cavity (1.5 • 10^^ ergsec^^Hz-^ for model BC at lOkeV). At 10 MeV 
the emitted power from the cavity is 2.7- 10^^ergsec~^Hz~^ and is dominated by the hotspot 
which contributes 8.0 ■ 10^"^ ergsec~^Hz~^. 

5.5 Caveats 

Which are the astrophysical and numerical limitations of our simulations? 

We have used a uniform external atmosphere, but extragalactic jets after having crossed a 
galactic halo will propagate into a much more diffuse intergalactic (or intra-cluster) medium. 
Both the pressure and the density decrease with the distance from the galactic nucleus. A 
density declining atmosphere will cause a widening of the jet (and of the cavity), and a 
substantial deceleration provided that the jets are not inertially confined. The effects of 
(continuous or abrupt) density changes in the external medium and its impact on the long 
term evolution of Newtonian jets has been investigated, e.g., by Hooda et al.(Hooda et 
al.l994). Such simulations involve additional, not very well determined model parameters 
and are more complex making them less attractive from the computational point of view, in 
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particular as our grid resolution (zones per beam radius) is much larger than that in Hooda 
et al.(1994). 

Our simulations are purely hydrodynamic, i.e., they do not consider the effects of mag- 
netic fields. Strong magnetic fields might confine the shocked jet material in an extended 
nose cone preventing it to be continuously deposited into the cocoon. An episodic release of 
thermal material in the cocoon would have important consequences for the stability of the 
jet. However, it is quite unlikely that at kiloparsec scales strong magnetic fields do exist. 
Instead, jets are expected to transport weak and mainly randomly oriented magnetic fields 
(e.g., Ferrari 1998). 

The grid resolution employed (six zones per beam radius at the injection nozzle) is far 
from being appropriate to account for phenomena like mass entrainment from the external 
medium. However, it is sufficient to describe the gross morphological features and to cap- 
ture the average cocoon dynamics. One should consider that, in order to compute up to a 
maximum physical time, tend- the required run time is proportional to (cells/Rb)^. Hence, 
increasing the grid resolution by a factor of two requires about 1200 hours of computing 
time to reach tend — 6 • 10^ y, and 9600 hours to simulate a typical lifetime of a FRH source 
{tend — 10^ y). In order to check whether the chosen resolution was sufficient to capture, at 
least the qualitative morphology, we performed a set of shorter simulations tend — 2.3- 10^ y at 
different grid resolutions (Fig. 12). With four zones/i?6 the short acceleration phases caused 
by the vortex shedding are missing, while they are present with six zones/ Rb. At even higher 
resolution the qualitative changes are small. 

Finally, we have restricted ourselves to axisymmetric jet flows. This is a considerable 
restriction, because we force the jet to propagate along the symmetry axis. As, e.g.,Aloy 
et al.(2000) showed, three-dimensional jets may wobble and their effective head area may 
change appreciably. Consequently, two possibilities arise: (i) the effective area grows (as 
the result of, e.g., the dentist drill effect, Cox, Gull & Scheucr (1991)) leading to a faster 
deceleration , or (ii) the wobbling of the Mach disk reduces the cross sectional surface of 
the jet head leading to an acceleration (Aloy et al.2000). It is difficult to extrapolate the 
results of Aloy et al.(2000) and to decide whether or not such mechanisms operate in the 
long run, and what their consequences are. Thus, three-dimensional numerical simulations 
are required to clarify this issue. 

Because of the assumed axisymmetry we encountered numerical problems when sim- 
ulating under-pressured jets (see Sect. 2.2). These problems would not appear in three- 
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dimensional simulations as the jet fluid will bypass the blobs of external medium material 
trying to block the inlet. 

6 SUMMARY AND CONCLUSIONS 

We have performed the longest and best resolved large-scale simulations of relativistic jets 
up to date. The simulations extend up to an evolutionary age of 6 • 10^ y which is only a 
few times smaller than the ages of typical powerful radio sources. The simulated jets differ 
in their density, temperature and composition, but have similar gross properties (kinetic 
luminosity and thrust). Hence, they may represent a relatively early stage of the observed 
extragalactic jets and their radio lobes. We have considered three models including a hot 
and a cold jet made of a pure (leptonic) electron-positron plasma (models LH and LC) and 
one baryonic model BC consisting of an electron-proton plasma. The leptonic hot model LH 
has a factor of 100 lower density (77 = 10~^) than the other two models (77 = 10~^). 

Although there is a difference of about three orders of magnitude in the temperatures 
of the cavities inflated by the simulated jets, we find that both the morphology and the 
dynamic behaviour are almost independent on the assumed composition of the jets. Only 
the leptonic hot model LH behaves differently because of its very light beam. This reflects 
the well known influence of the density parameter rj according to which light jets are less 
stable than heavy ones. Model LH suffers from considerable mass entrainment into the beam, 
which eventually leads to the disruption of the beam. The instability grows on a time scale 
of a few million years, i.e., the evolution of model LH is inconsistent with the sizes and ages 
of the jets of typical powerful radio sources. In addition, the initial conditions chosen for 
model LH are quite extreme, and most likely cannot be found in actual jets. 

In all models two evolutionary epochs can be distinguished. During the first epoch multi- 
dimensional effects are not yet important and the jet propagates with a velocity close to the 
ID estimate. The Mach disk is only shghtly disturbed, the shedding of vortices is negligible, 
and both the hotspot pressure and the jet velocity are nearly constant. The aspect ratios of 
the cocoon and of the cavity increase with time, while the average cavity pressure decreases. 
The second epoch starts when the first larger vortices are produced near the jet head. The 
beam cross section increases and the jet decelerates progressively. The deceleration is re- 
peatedly interrupted by short acceleration phases caused by the temporary replacement of 
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the Mach disk by a conical shock. The velocity increase during the acceleration phases can 
exceed 100%. 

The evolution of the cocoon and cavity aspect ratios of all models is in agreement with 
an extended version of the simple theoretical model proposed by Begelmann & Cioffi (BC89) 
for the inflation of the cocoon. The cavity aspect ratio of the simulated models is too large 
according to current observations, i.e., the simulated jets are too elongated. This situation 
may change when the simulations are either extended even further in time or when the 
assumption of axisymmetry is relaxed. However, both approaches would require prohibitive 
amounts of computational resources and axisymmetry is, in general, a good approximation 
to study extragalactic jets. 

The size and shape of the cocoon depends on 77, i.e., on the ratio of the densities in 
the jet and in the ambient medium. According to our definition of the cocoon - the region 
containing material with a beam particle fraction 0.1 < X < 0.9 - wc find that a light 
jet (model LH) has a very thin cocoon restricted to a narrow layer girding the beam. Jets 
with heavier beams (models BC and LC) possess more extended cocoons, more similar to 
those of observed jets. The differences in the cocoon sizes become particularly evident when 
comparing the synthetic radio maps that we have computed for all models. The cocoon 
aspect ratio found for the leptonic hot model LH is unrealistically large, which provides 
another argument against its realization in nature (see above). 

The propagation velocities and the hotspot pressures of our jets agree very well with the 
analytical and observed values. The strength of the equipartition magnetic field comfortably 
matches the expectations for a young powerful radio source (whose magnetic fields are 
thought to be stronger than those of evolved radio sources). The beam velocities of the 
models are relativistic (F ~ 4) at kiloparsec scales. Provided this is not an artefact caused 
by the assumed axisymmetry, the results strongly support the ideas of Tavecchio et al.(2000) 
and Celotti et al.(2001) that the X-ray emission of several extragalactic jets may be due to 
relativistically boosted CMB photons. 

The radio emission of all models is dominated by the contribution of the hotspots. The 
two cold models (BC and LC) have very similar radio-morphologies with a continuous beam 
connecting the inlet with the hotspot. This is different from model LH where the synthetic 
radio maps show a discontinuous and knotty beam even at a viewing angle of 45° (where 
the effects of Doppler boosting are still small). 

We find that the thermal bremsstrahlung emission is strongly suppressed at very high 
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energies (~ lOMeV) in the leptonic cold model because the average temperature of its 
cavity (~ 10^'' K) is much lower than that of the two other models (~ 10^^ K). The baryonic 
and the leptonic hot model both show a gradual transition in the characteristics of their 
thermal bremsstrahlung emission. At low energies the emission is dominated by the shell, 
while the interior of the cavity dominates at high energies. All models exhibit a depression 
in the X-rays surface brightness of the cavity interior in agreement with recent observations 
(McNamara et al.2000). These observations also show a lack of strong X-ray emission from 
the boundaries of the radio lobes of Hydra A. This does not contradict our results as Hydra A 
is a FRI source, while we have been simulating powerful radio sources. In a powerful source 
the bow shock is a strong shock and, hence, it may produce substantial X-ray emission. FRI 
sources flare into the external medium and, most probably, they do not have strong bow 
shocks. 
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Figure 1. Values of rj (solid lines) and T;, (dashed lines) for leptonic jets (X;;, = 1) as a function of Rf, and Wi, after fixing 
Lkim f j'') ^ib ^'^^ the external medium according to Tab. 1. There are no physical solutions in the areas without Tj-contour 

lines. The thick, solid lines correspond to values of r/ (from top to bottom) of 0, 10^^, 10^''', 10^^. The thick, dashed lines 
correspond to values of Tj, (from left to right) of 10^"K, 10"K and 10*^K. For baryonic jets the diagram would be similar, 
except for the values at the thick, dashed lines which would be three orders of magnitude higher: lO-'^'^K, and lO^^K, 

respectively - the reason being that for the same value of i), one would have a value of Xn, three ten folds larger. The squares 
mark the situation of the three simulated models in this diagram (see 2.3). 



Figure 2. Snapshots of the logarithm of the density (logm (p/pm)) for the models BC (top panel), LC (central panel) and LH 
(bottom panel) at t ^ 6.3 ■ lO^'y. The black lines arc iso-contours of the beam mass fraction with X = 0.1 (outermost) and 
X = 0.9 (innermost). These values correspond to the boundaries of the cocoon and the beam, respectively. Please, note that 
the colored version of the figure is only provided in the electronic version of the journal. 
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Figure 3. Snapshots of the logarithm of the temperature (logi()(T), T in Kelvin) for the models BC (top panel), LC (central 
panel) and LH (bottom panel) at t Ri 6.3 • lO^y. The black lines are iso-contours of the beam mass fraction with values 0.1 
(outermost) and 0.9 (innermost), that correspond to the limits of the cocoon and the beam, respectively. Please, note that the 
colored version of the figure is only provided in the electronic version of the journal. 
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Figure 4. Snapshots of the Lorentz fa<;tor W for the models BC (left panel), LC (central panel) and LH (right panel) at 
t fs 6.6 • lO^y for model LH and t i« 6.3 • lO^y for models LC and BC. The black lines are iso-contours of the beam mass fraction 
with values 0.1 (outermost) and 0.9 (innermost), that correspond to the limits of the cocoon and the beam, respectively. Please, 
note that the colored version of the figure is only provided in the electronic version of the journal. 



Figure 5. Profiles along the jet axis (r = 0) at t ^ 6.3 ■ lO^'y of (a) density (in units of pm), (b) pressure (in units of pmC^), (c) 
specific internal energy density (in units of c^), (d) local sound speed (in units of c), (e) temperature (in Kelvins), (f) Lorentz 
fa<;tor and (g) beam mass fraction. The units given at the abscissas are in kpc. 



Figure 6. Time evolution of (a) the jet length Ij, average cavity radius Tc, average eoeoon radius rj, (b) the aspect ratios 
(length/radius) of eoeoon and cavity, (e) the hotspot pressure Pfi^, average cavity pressure Pc and (d) jet velocity Vj for models 
BC (dashed), LC (solid) and LH (dotted), respectively. The jet velocities of models LC and LH have been divided by 10 and 
100, respectively, in order to enhance the readability of the figure. 



Figure 7. Bow shock velocity Vc for the three models. 



Figure 8. Hotspot pressure P/j^ (solid) and ram pressure pm^c °f the bow shock (dashed). 



Figure 9. Ratio of beam particles to external particles in the cocoon. 

Figure 10. Simulated total intensity (in logarithmic arbitrary units) radio maps for models BC, LC and BH (from top to 
bottom) after 6 • 10® y, and for a spectral index of —0.7. Both, the jet (to the right) and counter-jet emission are shown for a 
viewing angle of 45° in all three models. To allow for a better comparison with observed sources the images have been convolved 
with a circular Gaussian beam with FWHM of lORf,. The same limiting "noise" level is used in all frames to allow for a better 
comparison of the models. 

Figure 11. Surface brightness per unit of frequency for model BC at lOkeV (top) and lOMeV (bottom), respectively. No plots 
are shown for the other two models, as the plots for model LH and the plot of model LC at lOkeV would be indistinguishable 
from those of model BC, and as model LC shows no bremsstralilung emission, at 10 MeV. The colour scales are normalised to 
maximum values of 4.5 • 10~^^ ergsec~^Hz~i cm~^ (top) and 2.8 • 10~^^ ergsec~^Hz~i cm~^ (bottom), respectively. Please, 
note that the colored version of the figure is only provided in the electronic version of the journal. 



Figure 12. Beam velocity averaged over cross sectional cuts through the beam perpendicular to the jet axis for different 
grid resolutions. The solid line corresponds to a resolution of 6 zones per beam radius (working resolution). The dashed and 
dashed-dotted lines correspond to test resolutions of 4 and 10 zones per beam radius, respectively. 
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